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Abstract: Charged-Current Quasi-Elastic (CCQE) neutrino scattering is the signal chan¬ 
nel for sterile neutrino oscillation experiments. Recent cross-section measurements have 
made it clear that the current understanding of this channel in the few-GeV region is in¬ 
complete, and several sophisticated theoretical models have been proposed to tackle this 
issue, although it is not clear which model best describes the global dataset. In this pa¬ 
per we argue that the current uncertainty surrounding CCQE cross-sections is a serious 
problem for experiments seeking to produce sterile neutrino limits. We perform a sterile 
neutrino analysis with published MINERuA data as an illustrative example. We highlight 
the need for caution in interpreting sterile neutrino limits given the context of incomplete 
cross-section model information. 



1 Introduction 


Accelerator neutrino experiments in the few-GeV region, with detectors at short baselines, 
are used both to constrain sterile neutrino mixing models and to measure neutrino-nucleus 
scattering cross-sections. As the measured quantities are event rates - the flux multiplied 
by the cross-section - the measurement of either relies on some assumption about the other. 

For a long time, relativistic Fermi gas (RFG) models [1] have been used to describe 
charged-current neutrino-nucleus scattering in generators [2], The only free parameter 
unconstrained by electron scattering data in these models is the axial mass, Ma, which 
was well-constrained to be Ma = 1.014 ± 0.014 GeV [3] from deuterium scattering [4] 
and pion electro-production data [5]. It was therefore believed that the Charged-Current 
Quasi-Elastic (CCQE) cross-section was well-understood; however, recent neutrino-nucleus 
scattering data on heavy nuclear targets have produced much higher cross-sections, and 
much higher axial-mass values in this simple cross-section parametrisation [6-10]. This 
discrepancy is thought to result from additional nuclear effects which are not included in 
the RFG models [11], This has led in recent years to the development of more sophisticated 
models to explain the incompatibility between datasets. These models differ significantly 
in their prediction of outgoing particle kinematic distributions, and as such, the state of 
neutrino-nucleus scattering cross-sections in the few-GeV region cannot be said to be well 
understood. 

In this analysis we investigate the effect that different cross-section models of the CCQE 
interaction channel have on the limits produced by a short-baseline muon-neutrino disap¬ 
pearance analysis using a 3+1 mixing model. The cross-section models investigated are a 
small range of those currently available in generators. The MINERVA CCQE cross-section 
data in neutrino and antineutrino modes [12, 13] is used as an illustrative example, though 
the conclusions of this work apply to any sterile neutrino measurements made with accel¬ 
erator neutrino beams in the few-GeV region. We show that the choice of cross-section 
model has a significant impact on the sterile neutrino confidence limits produced, and ar¬ 
gue that the current uncertainty on the CCQE cross-section makes sterile neutrino limits 
in this energy range difficult to interpret. This work builds on work done in [14, 15] to 
show that modifications to Ma in the RFG model can affect the neutrino limits produced 
by sterile analyses. It complements other work investigating the effect that uncertainties 
in the cross-section models have on the reconstructed energy [16], fitted limits on 5qp [17], 
and atmospheric mixing limits [18] in a three neutrino framework. 

Whilst a fake data study would have been equally valid for the purpose of this analysis, 
we chose to use public MINERVA CCQE cross-section data, as a sterile neutrino fit of this 
kind has not yet been performed on these datasets. The NuWro Monte Carlo event generator 
[19] was used to produce differential cross-sections from initial event rate predictions for 
the CCQE cross-section models detailed in Section 2. Sterile neutrino induced biases to 
these predictions were produced by folding in a muon-neutrino survival probability under 
the 3+1 mixing model described in Section 3. Each sterile hypothesis was then fitted to 
the MINERi'A dataset detailed in Section 4 and a x 2 statistic was calculated as described 
in Section 5. For each of the cross-section models we investigated, y 2 scans were performed 
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in the sin 2 2 6^ — Am^ plane. The resulting confidence intervals are discussed in Section 

6 . 

2 Cross-section models 

This section describes the key features of the cross-section models considered in this analysis. 
There are three nuclear models, described in Section 2.1, one of which is the familiar Smith- 
Moniz RFG model [1] used in many generators and past analyses. Models of additional 
nuclear effects are described in Section 2.2. 

2.1 Underlying nuclear model 

Dipole axial form factors [20] and BBBA05 modifications [21] to vector form factors were 
used consistently for all of the models described in this section. 

Relativistic Fermi Gas (RFG): Nucleons are treated as quasi-free with a nucleus- 
dependent Fermi momentum and constant binding energy, Ef, [1]. This model uses the 
impulse approximation where the neutrino interacts with one nucleon only. In the RFG 
model all states up to the Fermi momentum are filled, so interactions where the outgoing 
nucleon is not outside the RFG distribution are Pauli blocked. The Bodek-Ritchie modifi¬ 
cation to the RFG model is included, which adds a higher momentum contribution due to 
short-range correlations between nucleons [ 22 ], 

Benhar Spectral Function (SF): A nucleus-dependent description of nucleon kinematics 
within the nucleus, in terms of its removal energy and momentum [23]. Approximately 20% 
of the cross-section is due to short-range correlations of nucleons (quasi-deuterons). The 
impulse approximation is used consistently; the interaction is with a single nucleon even for 
correlated states. Pauli blocking is approximated by a nucleon-dependent cut-off [24], 
Local Fermi Gas (LFG): Similar to the RFG but the binding energy, £ 5 , varies with the 
nucleon position within the nucleus, producing a more realistic Pauli blocking effect [25, 26]. 


2.2 Nuclear effects 

Recent models attempt to explain the large MiniBooNE axial mass value in terms of mod¬ 
ifications to CCQE interactions within the nucleus. We consider three such enhancements 
to the standard model. 

Transverse Enhancement Model (TEM): A four-momentum transfer dependent mod¬ 
ification to the magnetic form factor [27]. The modification is obtained by fitting to an 
experimentally observed excess in the ratio of transverse to longitudinal quasi-elastic re¬ 
sponse functions from electron scattering data [28]. 

Random Phase Approximation (RPA): A modification to the quasi-elastic propaga¬ 
tor, which accounts for long range nucleon-nucleon correlations within the nuclear medium 

[291 ■ 

Nieves multi-nucleon interaction model: A microscopic model that sums over possible 
W boson absorption modes, where the interaction is with two or three nucleons [30]. Note 
that this is an explicit contribution beyond the impulse approximation which has final states 
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that are largely indistinguishable from CCQE interactions (CCQE-like), and will therefore 
enhance reported CCQE cross-section measurements. 

2.3 The Nieves model 

The Nieves model [30] is a consistent description of the CCQE-like cross-section which 
incorporates the LFG, the RPA, and the Nieves multi-nucleon interaction model. From 
now on the “Nieves model” will refer to this combination. 


3 Sterile models 


3+1 neutrino models extend the 3x3 PMNS matrix by including an additional, predom¬ 
inantly sterile, mass state which is heavier than the other three neutrinos. Over short 
baselines, the three active mass states can be approximated as degenerate, and a two neu¬ 
trino mixing equation can be used to describe mixing between any of the three active states 
and the larger sterile mass state. The survival probability of a muon (anti-)neutrino can 
then be calculated using [31]. 


P (V. -VV ) = 1 - sin^ 26 


MM 


sm 


1.265 Am| 4 [eV 2 ] L[km] 
E v \ GeV] 


(3.1) 


where 


sin 2 20w = 4(1 - |[/ m | 2 )IM 2 ■ (3.2) 

4 MINERVA CCQE data 

This analysis uses the MINERVA v p and CCQE cross-section measurements [12, 13]. 
The data were taken on a CH target and are presented as a differential in reconstructed 
four-momentum transfer, Qqe- The key experimental details are summarised in Table 
1 [32], The public data release includes the full covariance matrix including correlations 
between the two datasets. 

The reconstructed neutrino energy, E„ , and four-momentum transfer, Qqe, are de¬ 
rived from the outgoing lepton kinematics (E^^p^, 6 p ) and the measured target binding 
energy E & for the target by assuming the pure two-body kinematics of the RFG model: 

e q E = m n ~ ( m p ~ E b) 2 ~ml + 2 (m p - E b )E^ 

" 2(m p — E b — Efj, + cos 6 p ) 

Qqe = 2 E^iE,, - p p cos 6 P ) - m 2 . (4.2) 


5 Fitting method 

We used NuWro to make Monte-Carlo (MC) comparisons with the MINERVA datasets 
for each of the cross-section models. Sterile neutrino induced biases were introduced by 
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Neutrino Run 



Distance to target, L (km) 

1.04 

1.04 

Energy range (GeV) 

1.5 <E U < 10.0 

1.5 < E u < 10.0 

Protons on target (POT) 

1.014 x 10 20 

9.42 x 10 19 

Integrated flux (y cm -2 POT -1 ) 

2.429 x 10 -8 

2.916 x 10 -8 

Target material 

CH 

CH 

Binding energy (MeV) 

30 

34 


Table 1 . Specifications of the MINERVA datasets used in this analysis. The cross-correlations 
between the neutrino and antineutrino datasets provided in refs [12, 13] allowed both datasets to 
be fitted simultaneously. The distance to target was approximated as the distance from the NuMI 
target to the MINOS near detector [33, 34]. 


re-weighting the flux. This approach allows large samples to be generated with minimal 
computational overhead. A y 2 minimization using the MINUIT [35] fitting package was 
used to determine best fit sterile parameters and calculate limits in the sterile mixing plane. 

MC events for each cross-section model were initially generated with a flat true neutrino 
energy (E v ) distribution across the experimental range. 

The effect of a sterile neutrino is to modify the effective flux since sterile neutrinos do 
not interact. We re-weight the effective MINERVA flux according to the survival probability 
for a given sterile hypothesis and recalculate the derived cross-section according to steps 
1-6 shown below. 


1 . Events were binned into a histogram R(Qq E , E u ) where Qq E was calculated using 
Equation (4.2). 

2. R(Qq E , E v ) was normalised to the total MC cross-section a MC . 

3. R(Qq E , E v ) was weighted to the published MINERVA flux distribution 4>. 

4. The survival probability for the j th E u bin, P(Am 2 4 , d E u, Ej), was calculated by 
averaging Equation (3.1) over 40 equally spaced points within the bin. 

5. R(Qq E , E u ) was multiplied by P (Am| 4 , 0^, Ej) introducing a sterile bias. 

6 . A cross-section histogram B(Qq E ) was created by projecting R(Qq E , E u ) onto the 

Qqe axis - 

The i th Qq E bin of B(Qq E ) is thus given by 


Bi = da f = 

3 


R 




EfcZ R kl 


xa MC x 


®(Ek) X ^ (^ m 24> Ej) 


(5.1) 


To reduce the statistical error from the MC to negligible levels, a large number of events 
(10 7 ) were generated for each cross-section model. The statistical error on each Qq E bin is 
then less than 0.1%. 
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Figure 1. Neutrino (left) and antineutrino (right) flux-averaged cross-section predictions are shown 
for all cross-section models investigated without any sterile neutrino bias. The upper panel shows 
the differential cross-section and the lower panel the ratio of model to data. In the ratio model 
and data are area normalised. The error bars on the MINERVA data include both statistical and 
systematic errors. 

Unbiased cross-section predictions corresponding to the null hypothesis are shown for 
each cross-section model in Figure 1. The effect of sterile neutrino induced biases on the 
RFG + TEM model over a range of mixing parameters can be seen in Figure 2. Changes in 
E u distributions introduced by the bias have only a small effect on the shape because the 
peak neutrino energy is higher than the experimental kinematic limit 2 m p E u > Q 2 [36]. 
The shape’s response to sterile modifications is likely to be larger for other experiments. 

The initial y 2 definition used to fit the sterile hypotheses to data is 

16 16 

x’-EEO'"-". D ) V("f -T) ■ <«) 

i=l j=1 

where Mij is the covariance matrix, u[ J are the measured differential cross-section in Qqe 
bins, and are calculated using Equation (5.1). It was found that minimizing this statistic 
gave results far below the data points, an effect consistent with “Peelle’s Pertinent Puzzle” 
(PPP) [37], 

PPP can occur when fitting to a dataset containing large correlated uncertainties be¬ 
tween all bins. If the total normalisation is reduced in the fit, the relative size of the shape 
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Figure 2. Example sterile hypotheses for the RFG+TEM antineutrino model. Ratio to shape- 
only data (area normalised to unity) is shown to highlight the limited shape sensitivity to sterile 
modifications. 

errors increases, thus appearing to improve the agreement even if the shape of the predic¬ 
tion has not changed [37]. This causes the fit to prefer parameter values which predict a 
distribution that lies far below the data [38]. We avoid the PPP problem by separating the 
MINERVA covariance matrix into a total normalisation error, e = 10.9%, and a shape-only 
matrix, M'^ lape [39]. 

The extracted shape-only covariance matrix, M^' ape . could not be inverted analytically 
as a result of rounding errors in the published MINERz/A data. We dealt with this problem 
by using the two-step method of ref [40] , in which the bin errors and correlations are treated 
separately. The alternative x 2 definition is given by 



where a % are the uncorrelated shape-only statistical errors for the dataset. The shape- 
only correlated uncertainties are contained in Ci(A~ 1 )ijCj which is defined in Appendix 
A. The advantage of this procedure is that A t j is an invertible matrix. The constrained 
parameter a normalised the theoretical predictions to the total measured cross-section in the 
experimental range, allowing the square bracket to represent shape-only contributions while 
the final penalty term reflected the difference in normalisation between the MC and data. 
The combination of these techniques was found to be a robust way to fit a highly correlated 
dataset affected by PPP (for a more detailed explanation of PPP and the definition of (5.3), 
see Appendix A). 

For each cross-section model, parameter scans were performed for values in the ranges 
0.0 < sin 2 2 6up < 1.0 over 500 evenly spaced bins and 0.1 eV 2 /c 4 < Am| 4 < 100 eV 2 /c 4 
over 1000 logarithmic bins. For each parameter bin a sterile bias was introduced using the 
bin centre co-ordinates and a \' 2 value calculated using Equation (5.3). The minimum x 2 
values from the scans were passed as starting assumptions to M1NU1T which then found a 
true x 2 minimum in the parameter space [35, 41]. Best fit points and minimum y 2 values 
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Model 

RFG 1.35 

RFG 0.99 

TEM 0.99 

SF 0.99 

NEV 0.99 

Nucleon distribution 

RFG 

RFG 

RFG 

SF 

LFG 

M a (GeV/c 2 ) 

1.35 

0.99 

0.99 

0.99 

0.99 

Enhancements 

- 

- 

TEM 

- 

Nieves T RPA 

Null x 2 /15 

2.332 

2.433 

1.663 

2.833 

2.971 

Best sterile x 2 /13 

1.803 

2.803 

1.628 

3.253 

2.943 

Best sin 2 20 /t/J 

0.638 

0.817 

0.322 

0.000 

1.000 

Best A77124 (eV 2 ) 

8.463 

0.370 

5.913 

0.104 

1.073 


Table 2. Null hypothesis and best fit sterile hypothesis % 2 /Ndof values. 



Figure 3. 90% CL mixing parameter contours and best fit points (starred) for the cross-section 
models investigated. 

for each model can be found in Table 2. A Ax 2 method was used to produce 90% CL 
confidence limits around the best fit points as shown for all cross-section models in Figure 
3. The Id confidence limits for the separate shape-only or normalisation penalty terms 
from Equation (5.3) are compared in Figure 4 for the Nieves and SF models to highlight 
the relative strength of the normalisation term. 

6 Discussion and conclusion 

The X ' 2 values in Table 2 and the contours in figures 3 and 4 demonstrate the sensitivity 
of sterile neutrino fits to the adopted cross-section model. Some models, e.g. RFG 1.35, 
exclude the null hypothesis at >99%CL, while others, e.g. SF 0.99, prefer it. The choice 
of axial mass value is a critical parameter (compare RFG 1.35 with RFG 0.99), but other 
features of the models tested also have significant effects (consider RFG 0.99 and TEM 
0.99). It is clear that in many situations the choice of cross-section model can completely 
dominate the results obtained in sterile neutrino searches. 

In the case of the MINERVA dataset considered in this study, the weakness of the 
final limits can be attributed to the large normalisation error. It is worth noting that the 
magnitude of the disagreement between models is likely to increase when analysing data 
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Figure 4. Overlaid shape-only, normalisation-only, and combined contours at the la confidence 
limit. The Nieves (left) and SF (right) models are shown. Best fit points are indicated by a star. 

with a smaller normalisation error or a stronger shape response to sterile neutrino biases. 
We conclude that sterile mixing limits obtained in this way are subject to large systematic 
uncertainties, until they can be repeated with a well-motivated theoretical model that agrees 
well with existing neutrino data. 

A x 2 definition and fake data study 

There is a well documented problem that can arise when fitting data with a covariance 
matrix that contains large correlated uncertainties between bins, as in Equation (A.l) [39]. 
By suppressing the normalisation of the prediction the y 2 is reduced, leading to a best fit 
distribution well below the data. This occurs because the covariance matrix is evaluated at 
a single point and the shape-only errors do not scale with normalisation. This problem is 
known as “Peelle’s Pertinent Puzzle” (PPP) [37]. 

x 2 = £ £ - v i) (^f - v f) ( AA ) 

i= 1 3=1 

Protecting against PPP is particularly important for sterile neutrino analyses, where a 
signal would involve a suppression of the overall normalisation, and a large normalisation 
uncertainty due to uncertainties in the flux prediction is common for accelerator exper¬ 
iments. PPP can be overcome by redefining the y 2 in terms of the shape-only matrix, 
Msh°pe. an( ] scaling the total integrated MC cross-section to match the total integrated 
cross-section in the data. This definition effectively stops the fit from inflating the relative 
size of the shape-only errors [39]. 

16 16 

x 2 = £ E b" --A ( M Ae) ,.hf-T) (a.2) 

i= 1 3=1 13 

We obtain M^ ape and the total normalisation uncertainty e from the published matrix 
Mij using the MiniBooNE matrix seperation method reproduced in Equation (A.3) [42]. 
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(A.3) 


The matrix separation method involves summations over many matrix elements which 
can lead to large rounding errors in M^ lape if is given to limited precision, as is the 
case for the MINERVA data release. This can cause problems when inverting M^ ape . We 
modify the x 2 definition in Equation (A. 2) to avoid inverting the matrix M^ ape using the 
method given in ref. [40]. The final y 2 definition used in our fits is given by 



(A.4) 


where cr,; is the uncorrelated shape-only statistical error on the i th bin, and is the 
correlated shape-only systematic uncertainty between the i th and k th bins which can be 
calculated using P2k=i ^ik^kj = M s P ape — af . The vector C and matrix A are defined 
by 


C ‘ = E 


AflkfKP - (l'"/o)) 


err. 


-A-ij — Sij + 


Aik^ 

k= 1 


jk 


ert 


(A.5) 


The matrix A, which is inverted in this method, is less susceptible to the rounding error 
problems than the full matrix MP ape . Note that with more complete information the 
matrix A can be reduced to a N x N matrix for N systematic errors as described in ref. 

[401- 

To test the y 2 statistic as a method of fitting for sterile induced biases, we conducted 
a fake data study using the RFG 1.35 dataset. We generated 30 sets of sterile parameters 
and generated fake data for each using the following method. 

First, MC was generated with the given parameter values according to the method 
described in Section 5. Then the systematic covariance matrix was calculated from the 
MINERVA matrix, M t j , and added to a diagonal matrix representing the statistical variance 
for the fake data study to create a fake data covariance matrix, M?, as shown in Equation 
(A.6). This covariance matrix would have been produced if the fake data reflected nature, 
and was measured in the MINERVA detector. 



Z'fVf 

vP v P 

* J 


[ Mij ~ S i3 ((O 2 - (a"°)?)] 


(A.6) 
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Throwing the matrix using the Cholesky decomposition method and adding the result to 
the nominal fake data produces a realistic fake data sample [43]. The residuals from 2000 
throws were calculated and fitted with a Gaussian to look for biases for each of the 30 sterile 
parameter sets investigated. These were found to have pulls away from the true parameter 
in the range of —0.05 to 0.10 for sin 2 2 9^ and —0.24 to 0.28 for Am| 4 , and widths in the 
range of 0.81 to 0.98 for sin 2 2 9^ and 0.75 to 0.98 for Am| 4 . It was concluded that the y 2 
statistic is a good estimator of the central value and was appropriate for the analysis. 
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